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Inviscid  stability  analysis  has  been  applied  to  the  mixing 
layer  profile  of  an  axisymmetric  jet  and  a coflowing  stream.  A 
collection  of  computer  subprograms  has  been  developed  to  solve 
the  resulting  eigenvalue  problem.  The  effect  of  changing  the 
velocity  profile  and  its  parameters  can  be  easily  assessed. 
Results  for  Gaussian  profiles  are  included. 
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Introduction 

The  mixing  region  between  a jet  and  a coflowing  stream  and  the  resulting 
large-scale  structures  are  currently  being  studied  at  NBS  using  a time- 
dependent  axisymmetric  computer  code.  When  these  structures  pair  and  merge, 
there  is  improved  mixing  between  two  streams,  which  is  useful  in  practical 
applications.  In  order  to  better  understand  the  flow  development,  controlled 
forcing  can  be  applied  and  the  effect  of  certain  types  of  perturbations  on 
the  growth  of  the  mixing  layer  observed. 

Linear  stability  theory  appears  useful  in  predicting  some  behavior  of 
the  large-scale  structures  in  a jet  [1].  The  evolution  of  vortices  seems 
to  be  related  to  the  unstable  frequency  range  determined  for  the  upstream 
velocity  profile.  Controlled  perturbations  applied  at  the  most  unstable 
frequency  and  its  subharmonics  can  affect  vortex  merging  and  the  growth  of 
the  mixing  region. 

In  the  present  work,  the  jet's  mixing  layer  is  studied  using  stability 
analysis.  The  objective  is  to  determine  the  spatial  stability  of  the  jet 
profile  and  use  the  results  as  the  controlled  perturbations  at  the  inlet  of 
the  jet  computation.  Although  the  procedure  developed  applies  for  arbitrary 
profiles,  the  velocity  ratio  of  the  jet  and  freestream  has  been  the  main 
parameter  varied.  The  analysis  assumes  that  the  flow  is  incompressible, 
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inviscid,  and  axi symmetric,  and  that  linearized  equations  are  applicable. 
The  validity  of  the  results  can  be  determined  by  comparisons  of  computa- 
tions and  experiments. 

In  the  following  section,  the  equations  used  in  the  stability 
analysis  will  be  presented.  The  solution  technique  using  software 
currently  available  on  the  NBS  UNIVAC  1100/82  will  then  be  described. 
Results  for  Gaussian  profiles  are  also  included. 


Basic  Equations 


The  method  presented  here  parallels  the  solution  procedure  described 
in  reference  [2].  The  axisymmetric  Euler  and  continuity  equations  for  an 
incompressible  fluid,  written  in  terms  of  the  three  variables,  u,  v,  p,  are 
the  basis  of  the  formulation,  u and  v are  velocity  components  in  the  axial 
(z)  and  radial  (r)  directions,  respectively,  in  cylindrical  coordinates, 
and  p is  the  ratio  of  pressure  to  constant  density.  Small  disturbances 
are  added  to  a known  solution:  u = U(r)  + u' , v = v' , p = p 1 , where  the 
primed  quantities  represent  the  disturbances.  Substituting  these  into 
the  equations  and  neglecting  higher  order  terms  involving  u'  and  v'  yield 
the  linearized  disturbance  equations.  Assuming  parallel  flow,  the  form 
of  the  disturbance  is  that  of  a traveling  wave  whose  amplitude  varies  with 
r and  which  moves  parallel  to  the  z-axis. 

- *(r)ei(aZ  - Bt) 


where  <j>  represents  u,  v,  or  p,  a is  the  wavenumber,  and  3,  the  frequency, 
of  the  disturbances.  Substituting  this  form  into  the  linearized  distur- 
bance equations  leads  to  three  ordinary  differential  equations  with 
complex  coefficients. 


M - e)v  = i ^ 
(all  - B)u  - i V 


ap 


For  a spatially  growing  disturbance,  a is  complex  (a  + ia- ) and  3 real. 
Eliminating  u and  v from  equations  (2)- (4)  gives  an  equation  for  the  pres- 
sure disturbance. 

d£ 
dr 


+ [-  - 
Lr 


2 dU 

(u  -de/a)]  - a2p  = 0 


The  boundary  conditions  to  be  satisfied  are  that  p remains  finite  at 
r = 0 and  p becomes  zero  as  r -*■  ».  Equation  (5)  and  the  boundary  condi- 
tions constitute  an  eigenvalue  problem.  The  objective  is  to  find  the 
corresponding  values  of  a and  s and  associated  eigenfunction  p(r)  for  a 
given  profile  U(r)  and  dll(r)/dr. 
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Computational  Procedure 


Several  computer  subprograms  are  used  in  finding  a solution.  The 
velocity  profile  being  used  in  the  present  study  is 

c - (r  - 1)» 

U = 1 - — — e 2gz  , r < 1 

/2lT 


(r  - 1)2 

u = LI  - (IL  + — - 1 ) e Zoz 
0 0 /27 


r > 1 


where  all  lengths  are  nondimensional ized  by  the  jet  radius  and  all  veloci- 
ties by  the  jet  centerline  velocity.  Equation  (6)  represents  either  a 
single  or  double  Gaussian,  depending  on  the  choice  of  U and  C//2tt".  U(r) 
and  dll(r)/dr  are  continuous  at  r = 1.  The  profile  U(r)°is  specified  in  the 
program  as  a function  UINLET(R)  and  dll/dr  as  DUDR(R).  Thus,  the  parameters 
or  form  of  equation  (6)  can  be  easily  changed. 

Equation  (5)  is  integrated  using  the  subroutine  CDRIV1  from  CMLIB 
(NBS  Core  Math  Library),  which  solves  initial  value  problems  for  complex- 
valued ordinary  differential  equations,  integrated  with  respect  to  a 
single,  real,  independent  variable.  Equation  (5)  is  expressed  as  a pair 
of  first-order  equations  in  subroutine  F. 


yi  = y2 


y2  = [■  f + 2 ^7  / (U  - 3/a)l  y2  + a2y i 


(6) 


(7) 


where  yl  = p and  y2  = dp/dr.  Subroutines  USERS,  JACOBN,  FA,  and  G,  which 
are  options  for  the  CDRIY  package,  are  included  as  dummy  routines  to  avoid 
error  messages  during  the  MAP  processing. 

Boundary  values  are  needed  at  r = 0 and  at  the  largest  value  of  r. 

At  each  of  these  locations,  dU/dr  vanishes  and  equation  (5)  becomes  a 
modified  Bessel  equation.  The  solutions  consistent  with  the  boundary 
conditions  are  p = a I (ar)  at  r = 0 and  p = b K (ar)  as  r becomes 
large.  This  leads  to  Boundary  values  at  r = 0 or  p = a 1(0)  = a, 
dp/dr  = act  I -.  ( 0 ) = 0.  Although  this  result  can  be  determined  directly 
from  axi symmetry,  the  Bessel  function  expression  is  used  to  obtain  the 
limiting  value  of  1/r  dp/dr  at  r = 0 needed  for  equation  (7).  Since  a 
is  arbitrary,  set  y^  = 1,  y2  = 0,  and  integrate  from  r = 0 to  r = 1 to 
get  the  inner  solution  f^(r).  The  integration  is  done  in  radial  incre- 
ments chosen  to  provide  velocity  values  at  the  grid  points  used  in  the 
jet  computation,  with  the  internal  step  of  the  integrator  adjusted  to 
just  reach  each  point. 
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The  initial  values  at  the  largest  value  of  r,  r , are  p = bK  (or  ), 
dp/dr  = - baK^ar  ).  Since  r is  large  and  evaluating  the  functions3* 
and  with  complexarguments  is  not  straightforward,  the  numerical  values0 
are  determined  from  asymptotic  expressions.  With  b also  being  arbitrary, 
set  yl  = K (ar  ),  y2  = - aMarmax)>  and  inte9rate  from  r = r__„  to 
r = 1 , as  aescrToed  above,  to  get  tne  outer  solution  f (r). 


max 


The  required  matching  conditions  at  r 
continuous.  This  reduces  to 


= 1 are  that  p and  dp/dr  are 


df.(D 

yu  1 


dr 


df0n) 

-Vu-ar- 


= 0. 


(8) 


The  expression  on  the  left  side  of  equation  (8)  is  written  as  a complex 
function  FF  whose  root,  a,  is  to  be  found.  The  integration  procedure  is 
included  in  the  function  subprogram. 

The  values  of  a are  determined  by  using  the  subroutine  ZANLYT  from 
IMSL  (International  Mathematical  and  Statistical  Libraries),  which  finds 
zeros  of  a single  nonlinear  complex  function.  For  each  value  of  6 chosen, 
ZANLYT  is  called  to  find  one  root,  with  one  guess  and  no  known  roots 
input.  The  first  guess  for  a found  most  appropriate  is  6(1  - i/2). 

When  the  root  is  returned,  several  tests  are  made  before  it  is  accepted. 
The  error  parameter,  IER,  must  be  zero,  indicating  convergence  was 
obtained.  For  the  problem  considered  here,  other  requirements  are  that 
a is  positive,  a.  is  negative,  and  the  phase  velocity,  c , = 6/a  , 
converges  to  one  it  6 = 0.  The  last  condition  is  checkedpGy  finding  a 
for  a slightly  smaller  value  of  6(a6  = 0.01),  using  the  first  a as 
a guess,  and  comparing  c , values.  If  c , is  less  than  one  and  decreases, 
or  c . is  greater  than  one  but  less  thanpT.5  and  increases,  the  value  of 
6 isprurther  decreased  and  a determined,  to  see  if  the  trend  reverses  in 
a short  interval  (20  A6,  A6  = 0.04).  If  the  first  root  is  rejected,  a 
second  guess  is  tried.  When  convergence  was  obtained  on  the  first  guess 
with  a positive  but  a.  not  negative,  the  complex  conjugate  of  a is 
guessed  for  the  second1 trial . Otherwise  the  guess  used  is  6 + 1/4 
+ i / 2 ( 1 - 6)  for  6 > 1 or  6(1  - i/10)  for  6 < 1.  For  the  cases  tested, 
a root  satisfying  the  conditions  described  above  was  found  without- 
additional  trials  needed. 

When  determining  a as  a function  of  6,  the  testing  is  done  only  for 
the  initial  value  of  6,  with  A6  = 0.04  in  the  c , comparison.  The  guess 
for  successive  points  is  the  value  of  a from  thipreceding  6 value.  This 
case  is  of  interest  for  finding  the  most  unstable  frequency,  the  value  of 
6 for  which  a-  has  the  largest  negative  value.  The  curve  is  started  from 
6=4  and  proceeds  to  smaller  values  of  6,  with  A6  = 0.04.  If  the  value 
of  a-  is  decreasing  as  6 decreases,  the  curve  is  continued  in  the  same 
direction  until  the  minimum  value  of  a-  has  been  passed.  If  increases 
in  the  first  set  of  intervals,  the  curve  is  continued  similarly  in  posi- 
tive increments  from  6=4  until  the  minimum  is  located. 
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Once  a has  been  found  for  a particular  value  of  g,  the  corresponding 
velocity  eigenfunctions  can  be  determined.  The  outer  pressure  solution 
f (r)  and  df  (r)/dr  are  multiplied  by  the  ratio  of  the  constants  b/a  = 
f°(l)/f  (1)  8 r df. (l)/dr  / df  (l)/dr,  derived  from  the  matching  conditions. 
Or*ie  arbitrary  constant  remains,  multiplying  the  entire  solution.  The 
velocity  v can  be  computed  at  each  value  of  r where  the  integration  was 
done  from  the  values  of  dp/dr  and  equation  (2).  The  u velocity  is  then 
calculated  from  equation  (3)  using  v and  the  p solution.  In  the  jet 

computation,  v is  needed  at  points  halfway  between  those  used  in  the 

integration,  so  averaging  is  done  there.  The  symmetry  conditions  at 
r = 0 give  u(-  Ar)  = u(Ar),  v(0)  = 0,  v(-  Ar)  = - v(Ar). 

From  equation  (1),  the  real  velocity  disturbance  at  z = 0 is  u‘  = 
u (r)  cos  3t  + u.(r)  sin  3t,  where  u = u + iu. . At  any  value  of  r,  u1 

i£  a maximum  when  Bt  = tan"1  (u./u  ),  giving  u'max  = /u|  + u.2.  The 

u and  u.  arrays  are  searched  to  determine  theradial  locations  where 
each  has  the  largest  magnitude.  At  these  two  values  of  r,  the  quantity 
0.1  U(r)/u*  is  computed  and  the  smaller  one  selected  for  the  constant 
a.  The  comprete  u and  v solutions  are  then  multiplied  by  this  factor. 

When  the  procedure  is  called  from  the  program  doing  the  jet  compu- 
tation, subroutine  PERTRB  is  used.  Inputs  to  this  subroutine  are  the 
array  of  r values  where  the  u velocities  are  needed  and  the  number  of 
points  in  the  array,  the  matching  location  (r  = 1),  and  S.  When  3 = 0 

is  input,  the  most  unstable  frequency  is  found  and  that  value  of  3 is 

returned.  Additional  outputs  are  a and  arrays  of  u , u . , v , and  v • . 

The  perturbations  applied  to  the  jet  velocities  are: 

N 

u'  = z a.  [u.  (r)  cos  ( 3^t  - a.  z)  + u.  (r)  sin  ( 3^t  - 

k=l  K r K r Ki  K 

N 

v'  = z a.  [v.  (r)  cos  (3tt  - a.  z)  + v.  (r)  sin  ( 3^t  - 

k=l  K Kr  K Kr  Ki  K 

where  N frequencies  have  been  included. 


Results 

Results  will  now  be  presented  for  the  inlet  velocity  profile  of 

equation  (6)  for  three  choices  of  the  constants:  (a)  a single  Gaussian, 

U = 0.65,  C//27  = 0.35,  (b)  a double  Gaussian,  U = 0.3,  C I JlH  = 0.8, 

and  (c)  a single  Gaussian,  U = 0.15,  C//2tt  = 0.85,  where  a = 0.1  in 

each  case.  These  profiles  are  shown  in  figure  1.  For  each  of  these, 

a was  determined  over  a range  of  3 * s . For  these  runs,  the  stopping 

criteria  in  ZANLYT  were  (1)  the  magnitude  of  the  function  FF  reaching 

or  falling  below  10"6,  (2)  two  successive  approximations  agreeing  in 

the  first  six  digits,  or  (3)  the  number  of  iterations  exceeding  the 

maximum  of  50.  The  most  iterations  actually  used  in  finding  any  one 

root  was  38.  The  requested  relative  accuracy  in  CDRIV1  was  set  at 

10"6.  The  value  of  r was  5.2. 

max 


°k z)^ 
r 

ak  Z)], 
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In  figure  2,  -a^  is  plotted  against  B>  where  the  variation  of  the 
peak  value  and  most  unstable  frequency  can  be  seen.  The  latter  ranges 
from  3.2  to  4.32  for  these  cases.  In  figure  3,  c , is  plotted  against 
B.  The  trend  toward  c , = 1 at  B = 0 is  evident  nere.  The  velocity 
ratio  listed  in  figuresn2  and  3 is  defined  as  1/U  . 


Concl us  ion 

A collection  of  subroutines  has  been  developed  to  solve  the 
eigenvalue  problem  resulting  from  applying  stability  analysis  to  the 
axi symmetric  continuity  and  Euler  equations.  The  velocity  eigen- 
functions are  directly  available  for  use  as  perturbations  to  the 
inlet  velocity  profile  in  the  jet  computation. 
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Fig.  1.  Velocity  profiles 
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